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Abstract 

This paper presents an algorithm to compute the Markov parameters of a backward- 
time observer for a backward-time model from experimental input and output data. The 
backward-time observer Markov parameters are decomposed to obtain the backward-time 
system Markov parameters (backward-time pulse response samples) for the backward-time 
system identification. The identified backward-time system Markov parameters are used 
in the Eigensystem Realization Algorithm to identify a backward-time state-space model, 
which can be easily converted to the usual forward-time representation. If one reverses time 
in the model to be identified, what were damped true system modes become modes with 
negative damping, growing as the reversed time increases. On the other hand, the noise 
modes in the identification still maintain the property that they are stable. The shift from 
positive damping to negative damping of the true system modes allows one to distinguish 
these modes from noise modes. Experimental results are given to illustrate when and to 
what extent this concept works. 


Introduction 

Identification is the process of developing or improving a mathematical model of a phys- 
ical system using experimental data to describe the input, output and noise relationship. A 
tremendous number of system identification techniques are available. Techniques to identify 
a model from input and output data typically contain two steps. First, a family of candidate 
models is chosen and then the particular member in this family is determined which satis- 
factorily describes the observed data based on some error criterion such as minimizing the 
measurement residuals due to the input and output noises. One of the most popular models 
is ARX model 1 where AR refers to the autoregressive part (related to output data) and X 
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refers to the exogeneous part (related to input data). In particular, the forward-time ARX 
model in the time domain is nothing but a linear finite difference equation describing that 
the current output measurement can be estimated from past output and input data. On the 
other hand, the backward-time ARX model means that the current output measurement can 
be estimated by the future output and input data. The ARX coefficients in the linear finite 
difference equation can be easily computed from input and output data. 

Insights into physical mechanisms of the system can usually be incorporated more easily 
into state-space models than the ARX models just described above. In the state-space form 
the relationship between the input, output, and noise data is written as a system of first- 
order differential equations using an auxiliary state vector. Note that it is difficult to directly 
compute a state-space model from input and output data. The state-space models can be 
constructed directly from the ARX coefficients using standard canonical forms including 
companions form, controllable/observable canonical form, Jordan canonical form, etc. One 
common problem is that an estimate of system order is required to be known a priori. 
Since the model order is generally unknown in practice, a trial and error process is required 
to find the best model order. In general, it results in a considerably overspecified model 
order particularly when the data contains considerable noise. Overspecification of the model 
order increases accuracy but discrimination of the system and computational portions of 
the models can be difficult. Recently, novel techniques referred to as backward methods 
using the backward-time model have been introduced in Refs. 2 and 3 to aid in the mode 
discrimination process. If one reverses time in the model to be identified, the damped true 
system modes become unstable, growing as the reversed time increases. On the other hand, 
the noise modes in the identification still maintain the property that they are stable. This 
is intuitively reasonable. If the data set is sufficiently long, an unstable noise mode would 
contribute to the pulse response data that grow unbounded as the time step in the data set 
increases. This is inconsistent with the expected contribution of noise in data. The shift from 
positive damping to negative damping of the true system modes allows one to distinguish 
these modes from noise modes. 

Aside from the problems of model order using standard canonical forms, numerical sensi- 
tivity is known to result from using the state space model whose elements are taken directly 
from the ARX model coefficients. The interested reader is directed to Refs. 4-6 for discus- 
sions of numerical sensitivities of companion forms, controllable/observable canonical forms, 
etc. The backward method introduced in this paper avoids numerically sensitive canonical 
forms by constructing a state-space realization from backward-time Markov parameters. This 
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method uses the observer-based approaches which are developed in Refs. 7-10 and referred to 
as the Observer/Kalman filter identification (OKID) method. The primary step is to use an 
observer to produce the backward-time Markov parameters (backward-time pulse response 
samples) from input and output data. The backward-time Markov parameters axe then used 
in the Eigensystem Realization Algorithm (ERA) 11-12 to construct a backward-time state- 
space model for its modal parameter estimation. As a by-product, a backward-time observer 
is obtained to be used for noise characterization. The identified backward-time state space 
model can be easily converted to the usual (forward-time) state-space representation. 

This paper starts by writing the relationship between the input and output histories in 
terms of the deterministic backward-time ARX model. Formulations are derived to compute 
the system backward-time Markov parameters and the backward-time observer gain from the 
backward-time ARX coefficients. The second part of this paper is to extend the deterministic 
case to the stochastic case with process and measurement noises. This will result in a 
backward-time ARMAX model where MA refers to the moving average part which is related 
to the residual of the measurement equation. The ARMAX model describes not only the 
deterministic but also the stochastic components of the response. The relationship between 
the backward-time observer and the Kalman filter is discussed. Experimental results in 
comparison with the forward OKID method 7-10 are given to illustrate the validity of the 
backward method presented in this paper. The experimental results are obtained from a 
10-bay truss structure having two accelerometers and two thrusters. 

Deterministic Backward-time ARX Model 

Consider a discrete multivariable linear system described by the state equation 

x(k + 1) = Ax(k) + Bu(k) (1) 


and the measurement equation 


y(k) = Cx(k) + Du(k) (2) 

where x(k) € Ft”, y(k) € R m , u(k) € R r , and k is the time index. Premultiplying Eq. (1) 
by A -1 and solving for x(k) yield 

x(k) = A _1 x(k + 1) — A~ l Bu(k) (3) 
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This is a backward-time model in the sense that the current state x (k) is determined by 
the future state x(k + 1) and the current input u(k). It is not a natual model and cannot 
be used in real time implementation. However, it can be very useful for modal parameter 
identification of a system model, which will be shown in the following. Substracting and 
adding the term Gy(k + 1) in Eq. (3) produces 


x{k) 


A~ l x(k + 1) - A~ l Bu(k) + Gy{k + 1) - Gy(k + 1) 

A~ l x(k T 1) — A~ 1 Bu(k ) + GCx(k + 1) + GDu(k + 1) — Gy(k + 


(A~ l + GC)x(k +!) + {- A' 1 B GD 


-G] 


u (k) 
u(k + 1) 

y(k + 1 ). 


1 ) 


(4) 


This is a backward-time observer model which observes the current state x(k) using the 
future state x(k + 1) and the observer gain G multiplied by a combination of inputs and 
outputs. The advantage of using the observer model instead of the original backward-time 
model is that the observer gain can be chosen so that the observer state matrix (,4 -1 + GC) 
is as stable as desired. Define 


A 

B 


A- 1 +GC, B= - A~ l B , 
[ B GD - G),v{k) = 


u (k) 
u(k + 1) 
y(k+l) 


Using the above definition, Eq. (4) can be shortened as 


(5) 


x(k) = Ax(k T 1) + Bv{k) 


( 6 ) 


Equation (6) looks identical to a typical linear discrete model except that it is backward 
in time. Prom Eqs. (2) and (6), the mapping between the input and output data can be 
established by the following equation 

y(k) = CA p x(k + p) + CA p ~ l Bv(k + p - 1) + ... + CBv(k) + Du(k ) (7) 


This equation computes an estimate of the current measurement y(k) using the p-step-future 
state x(k+p), current and future inputs, u(k ), . . . , u(k+p), and future output measurements, 
y(k + 1), y(k + 2), . . . , y(k -f p). Let G and p be chosen such that 

A i = (A- 1 + GC) 4 = 0; i > p (8) 


Physically, it means that the gain matrix G makes the backward-time observer state matrix 
A~ l + GC becomes deadbeat after p steps. It implies that the current output measurement 
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can be fully estimated by the future p inputs and outputs (see Eq. (7)) without any future 
state involved. The smaller the p is the shorter the number of the future inputs and outputs 
is needed to estimate the current measurement. However, there is a lower limit for the 
integer p which is related to the order of the system and the number of outputs. For the 
single-output observable case, there always exists a gain G which makes the state matrix 
A ~ 1 + GC becomes deadbeat as long as the chosen integer p is greater than or equal to the 
order of the system matrix A. For the multiple-output observable case, if p is chosen such 
that the product of p and the number of outputs is greater than or equal to the system order, 
a gain G always exists to make the observer state matrix A~ l + GC becomes deadbeat after 
p steps. The reader is directed to Ref. 10 for the proof of the above statements. 

Substitution of B and v defined in Eq. (5) into Eq. (7) yields 

y(k) + CGy(k + 1) + ... + CA p ~ x Gy(k + p) 

= (D + CB)u{k) + (i CAB + CGD)u(k + 1) -I- ... 

+{GA p ~ x B + CA p ~ 2 GD)u{k + p - 1) + CA p ~ x GDu(k + p) (9) 


or 

y(k) + any(k + 1) + ... + a p y{k + p) 



= Pou(k) + Piu(k + 1) + ... + Ppu(k + p) 

(10) 

where 

a i = CA i ~ l G-, i = l,2,...,p 



Pi = CA'B + CA'~ l GD\ i = l,2,...,p- 1 
#) = £> + CB; p p = CA p ~ x GD 

(11) 


Equation (10) is known to be a backward-time ARX model. At this moment, the relationship 
between the backward-time state space model, Eq. (3), and the backward-time ARX model 
becomes clear. Given the system matrices A,B,C,D, and the observer gain matrix G, the 
backward-time ARX mdel, Eq. (10), can be easily derived by using Eq. (11). On the other 
hand, if the coefficients o< and A (i = 1,2,..., p) of the ARX model are directly computed 
from the input and output data, are the system matrices A, B , C, D realizable from these 
coefficients? This question will be answered in the following. 
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Let Eq. (10) be rewritten in the following matrix form 


yik) - [ft, { -a p ) { dp— i -a p -i) 


{A -a,)] 


Now define 


V = 


u(k) 

f u(k + p) \ 

1 y(k + p)t 

(u(k + p- 1) 

\y(k + p-i) 

f u(k + 1 ) ) 
\y(k + 1) / 


} 


y = [»(0) y( l) 2 /( 2 ) 

Y = 


v(t - v ) } ; 


- 

a v } ( A>- i 

-a p _i} ••• { Pi 

«( 0) 

nil) 

uit-p) 

r u(p) \ 

r u(p+ 1) ) 

/«(01 

X y(v) J 

ly(p + i)J 

UWJ 

u(p - 1) 

1 / u(p) 1 

/«(*-!)! 

yip — i) 

J l yip) J 


/«(in 

/ «(2) \ 

(ui£-‘p+ 1) 

1 2/(1) J 

1 2/(2) J 

\yi£-p+l) 


}J 


( 12 ) 


(13) 


The matrix Y of dimension m x [r + (m + r)p] contains the unknown ARX coefficients ctj and 
Pi (i — 1,2, ... ,p), whereas the matrices y of dimension mx (t-p + 1) and V of dimension 
[r + (m + r)p] x (£-p + 1) contain the input and output data. Making use of Eq. (12) and 
the definition of Eq. (13) produces 

y=YV (14) 


which yields 


Y = y V' 


(15) 


where f means the pseudo-inverse of V. Equation (15) solves for the unknown ARX co- 
efficients £*j and Pi (i = 1 , 2 , . . . ,p). Having computed these ARX coefficients, the system 
matrices A, B, C, D can be constructed as follows. 

Using Eq. (11) yields 

D = -a~ l Ppi 

CB = po-D ; (16) 

CA { B = Pi-aiD 
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( 17 ) 


From Eqs. (11) and (16), one obtains the following Markov-like parameter sequences 

Y 0 =\CG CAG ••• CA^G] 

and 

Y a = [D CB CAB ••• CAp-'B] (18) 

All system information is embedded in the sequences Y a and Y a . The next step is to extract 
the system Markov parameters from the two sequences for the identification of the system 
matrices A,B,C,D. By simply applying the same approach as in Ref. 10, the sequences 
defined in Eqs. (17) and (18) can be decoupled to produce the following combined Markov 
parameters, 

Y C =[D C[B G] CA~ l [B G] ••• CA~ k [B GJ] (19) 

where k is an arbitrary integer which can be as large as desired. The sequence Y c is defined as 
the combined Markov parameters to describe the input-output map of a backward-time ob- 
server model as defined in Eqs. (4) in combination with the measurement equation, Eq. (3). 
The combined Markov parameter sequence can be used in conjunction with the ERA algo- 
rithm in Ref. 11 or 12 to realize a set of matrices 

[A~ l BC DG] 

Using the relationship B = — A~ X B , the forward-time representation can be computed 

[ABC DG) 

with B = —AB. The operator “inverse”, A~ l , shifts all the eigenvalues of A from inside the 
unit circle to outside the unit circle. In other words, the stable system modes will become 
unstable modes by the inverse operator. 

Stochastic Backward-time ARMAX Model 

Equation (11) is a deterministic backward-time ARX model representing the state space 
model given in Eqs. (1) and (2). In the following development, an extension to the case with 
process and measurement noises will be made. This will result in a backward-time ARMAX 
model that describes not only the deterministic but also the stochastic components of the 
response. Consider the system with process and measurement noises given below, 

x(k + 1) = Ax(k ) + Bu(k) + wi(k) 

y(k) = Cx(k) + Du(k) + tt /2 (k) ^ 
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where x(t) € R n , y(i) € R m , u(i ) € FT. Ideally, the process noise wi(k) and measurement 
noise w 2 (k) are two uncorrelated zero-mean stationary white noise processes. From Kalman 
filter theory, there exists a Kalman filter gain K such that the system given in Eq. (20) is 
equivalent to 

z(fc+l) = Ax(k) + Bu(k) + Ke(k) 

y(k ) = Cx(k) + Du{k) + e(k) < 21 ) 

where e(k) is a residual sequence which is white and uncorrelated with the measurements, 
and whose Euclidean norm is minimized. To derive a backward-time model for the above 
system, premultiply Eq. (21) by A~ x and solve for x(k), 

x(k) — A~ l x(k + 1) - A~ l Bu(k ) - A~ l Ke(k) (22) 

This is a backward-time stochastic state space model in the sense that the current state x(k) 
is described by the future state x(k + 1), the current input u(k), and the current residual 
e(k). If the original system is stable then the backward-time system will be unstable. To 
derive the backward-time ARMAX model, add and subtract the term Gy(k + 1) in Eq. (22), 

x(k) = A~ x x(k + 1) - A~ l Bu(k) - A~ l Ke(k) + Gy(k + 1) - Gy(k + 1) 

= A~ l x(k + 1) - A~ x Bu{k) - A~ l Ke(k) + GCx(k + 1) + GDu(k + 1) 

-\-Gs{k + 1) — Gy(k + 1) (23) 


With the definitions as given in Eq. (5), Eq. (23) can be simplified as 


x(k) = Ax(k + 1) + Bu(k ) + Ke(k) + GDu(k + 1) + Ge(k + 1) — Gy(k + 1) 

u(k) 


Ax(k+\) + [B GD -G] 


u(k -|- 1) 

lv(k + 1)J 


+ [K G] 


e(k) 

!.£(& + 1) J 


(24) 


where K = -A l K. The expression for the output y(k ) becomes 


y(k) = CAx(k + 1) + (CB + D) u(k ) + CGDu(k + 1) - CGy(k + 1) 

+ (/ + Ck) e{ k) + CGe(k + 1 ) (25) 


The estimated state vector x(k + 1) can be eliminated by successive substitution of Eq. (24) 
into Eq. (25). For example, 


£(fc + l) = Ax(k + 2) + Bu(k + 1) + Ke(k + 1) 

T GDu(k T 1) Ge(k T 2) — Gy{k + 2) 


( 26 ) 



Equation (25) then becomes 


y(k) = CA 2 x(k + 2)+{CB + D)u(k) + {cAB + CGD)u(k+l) + CAGDu(k + 2) 

- CGy(k+l)-CGy(k+l)-CAGy(k + 2 ) 

+ (/ + CK) e{k) + (CAR + CG ) e(k + 1) + CAGe{k + 2) (27) 

Repeating the substitution p — 1 times, and making use of the imposed deadbeat condition 
for G, i.e., 

A* = [A~ x + GG )' = 0 i > p (28) 

yields the following backward-time description relating the current and future outputs to 
current and future inputs and residuals, 

y(k) + CGy(k+l) + CAGy(k + 2)+ ••• + CA p ~ l Gy(k + p) 

= (D + CB) u(k ) + ( CAB + CGD ) u(k + 1) + 

+ (CA P ~ X B + CA P ~ 2 GD ) u(k + p - 1) + CA p ~ l GDu(k + p) + e{k) (29) 

where e(k) is a colored residual which is related to the white residual of the Kalman filter 
e(k) by the relation 

e{k) = (l + CR)e{k) + (CAR + CG)e(k+l)+ ••• 

+ (CA P ~ 1 R + CA p ~ 2 G) e(k + p- 1) + CA^Geik + p) (30) 

Define cti, 0i as in Eqs. (11) 

a* = CA'~ x G ; i = 1, 2, • • • , p 
fa = CA*- l B + CA'-'GD ; i = 1, 2, • •• , p-1 
0o = D + CB ; 0 P = CA p ~ l GD 

and 7 i as 

7< = CA'R + CA-'G ; i= 1, 2, p- 1 
7o = I + CR ; 'r P = CA p ~ 1 G 

Equations (29) and (30) become 

y(k) + atiy(k + 1) + ••• +a p y(k + p ) 

= 0ou(k) + 0iu(k + 1) -I- ••• + 0 p u(k + p) + e(k) (33) 


(31) 

(32) 


and 


e(k) = 7o£(fc) + 7i £ (k + 1) + • • • + 7i>£(k + P) 


(34) 
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Equation (33) above is now referred to as a backward-time ARMAX model. The coefficients 
of the backward-time ARMAX model are expressed in terms of the state space matrices 
A,B,C,D, a backward-time observer gain matrix G , and the (forward-time) Kalman filter 
gain K. Observe that the deterministic portion of the model described by the coefficients 
is the same as the backward-time ARX model derived in Eq. (10). 

System Identification using a Backward-time ARMAX Model 

The problem of system identification using this model structure consists of three basic 
steps: 

1. Determination of the model coefficients and 7 < from input and output data. 

2. Computation of the appropriate Markov parameters from the model coefficients 
and 7 j. 

3. Realization of the system state space model and the corresponding observer gains from 
the Markov parameters. 

Each of the three steps will be described below. First, it is convenient to rewrite Eq. (33) 
in the following form, 


y(k ) = [fa { 0p - a P } {0p- 1 -ap-i} 

e(k + p) 
e(k + p- 1) 

£(k + 1) 

For a given input-output time history, Eq. (35) can be written as 

y = YV + e 

e = AM/ + C (36) 

C = 7o£ 


+ [ip 7p-i ••• 7i 
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where y , Y, and V are given in Eq. (13), and e, £ e A, W are defined as 

e = (e(0) e(l) e(2) ••• e(£-p )] 

C = (C(0) C(l) C(2) ••• C (t-p)\ 

e = [e(0) e(l) e(2) ••• e(^-p)] 

A = [7 P 7 P -i ••• 7i ] 

e{p) £{p+ 1) • • • e(£) 

e{p - 1) e(p) • • • e(£ - 1) 

£(1) e(2) ••• e(t-p+ 1) 

At this point, the following observation can be made. The ordinary least-squares solution, 
which was derived from a purely deterministic consideration as given in Eq. (15), minimizes 
the norm of the colored residual e, whereas the true solution minimizes the white residual e. 
The least-squares solution, Eq. (15), is the solution of a linear problem, viz., y = YV. This 
is in contrast with the parameter estimation problem considered here, namely, 

y = YV + AW + 7o£ 

where in addition to the parameters in Y, neither A, 70 nor the residual in W are known. 
This is a non-linear problem since it involves products of the unknowns in KW and 7o£. When 
no closed-form solution is available, typical solution to this type of problem is iterative in 
nature. This approach has been considered previously in the forward-time observer/Kalman 
filter identification problem, which produces an improved solution under certain conditions 
when the ordinary least-squares solution is not optimal. The readers are referred to Eef. 
13 for further details. Examining the structure of Eq. (36) suggests the following iterat- 
ing procedure, which starts with the ordinary least-squares solution Yls obtained from the 
backward-time ARX approach. 

Step 1: Compute the least-squares solution Yls and the resultant least-squares error, 

= yv T {wry' (37) 

= y-Y LS V 

Step 2: Use the initial least-squares error eis as an initial estimate of the white residual 
e denoted by £(1), i.e., let £(q = e^s, and compute an estimate of backward-time ARMAX 
model coefficients in Y and A, denoted by Y(i), A(i), 

[V'oi, A(1)] = VmVfn {u m l fi))' 1 (38) 
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where 


Vd) = [ 2/(0) y{ 1 ) y{ 2 ) ••• y(£-2p)} 


U( 1) = 

' V ' 

£(i )(p) 

£ (i )(P + 1) 

£(i )(*-p) 

W(I) = 

£(i ){p- 1) 

£ (i)(p) 

••• £(1)(*-P-1) 


. £(i)(l) 

£(i)( 2 ) 

••• e ( i)(£-2p+l) 

The resultant least-squares error, denoted by £(i), is 

computed from 


C(D = J/( i) “ Y(i)V ~ A(i)W^(i) (39) 

Note that the parameters in Ym and A(i) do not include the first coefficient of the MA 
portion of the model denoted by 70. This coefficient is estimated in step 3. 

Step 3: Compute an estimate of the first coefficient of the MA portion of the backward-time 
ARM AX model, denoted by 

70 ,( 1 ) = C(1)4) [ € w € h)\ 1 ( 4 °) 

At this step, a complete set of the backward-time ARMAX model coefficient estimates is 
obtained in 70, (i), V(i), and A( ij. 

Step 4: Compute a new estimate of the whitened residual, denoted by £(2), from 

£ m - %,h)C(D ( 41 ) 

Step 5: Return to step 2 with the new whitened residual £(2) obtained in step 5 to compute 
new estimates of Y and A, denoted by Vj 2), A( 2 ), and then to step 3 and step 4 to compute 
new estimates of 70, denoted by 70,(2) , and a new whitened residual denoted by £(3). 

The above iterating cycle is repeated and the results examined after each cycle to de- 
termine a reasonable stopping point. Observe that after each cycle is completed, the data 
record to be used in the next cycle is shortened by p points. This is due to the backward-time 
nature of the ARMAX model which requires knowledge of p future data points to compute 
the current response. This is in contrast with the forward-time model where the initial 
condition can be neglected after the transient has decayed. In the backward-time model, 
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however, one cannot neglect future data points unless they are known to have decayed to 
negligible values. Specifically, given a data record of £ points, a backward-time model of 
order p can only estimate the first £ — p data points. Provided the data record to be used 
in the above iterating procedure is sufficiently long, the above shortening of the data record 
does not present a problem. 

One can compute the Markov parameter sequences of the backward-time observer and 
system 

Y 0 = [CG CAG ••• CAv-'G] 

Y s = [D CB CAB ••• CAv-'B] K ) 

by Eqs. (11) and (16). Furthermore, since 


CK =70 - / 
CA { k = 7i — a* 


(43) 


the following Markov parameter sequence can also be computed 


Y f = \CK CAk ■■■ CA p ~ l k] 


(44) 


The backward-time system and observer information is contained in the sequences Y a and Y s> 
whereas the noise information is contained in the sequence Yf. The rest of the identification 
procedure parallels the previous development with the exception that the identified noise 
model is being carried in the sequence Yf. Using the sequences defined in Eqs. (42) and 
(44), the following combined Markov parameter sequence can be computed, 

Y C =[D C[B G K) CA~ X B ••• CA~ k [B G K}\ (45) 


where A; is an arbitrary integer which can be as large as desired. The extra Markov parameters 
are computed by invoking the prescribed deadbeat condition 

CA { \B G K] = 0 i>p (46) 

Application of the Eigensystem Realization Algorithm (ERA) to the sequence defined in Eq. 
(45) will produce a realization of the backward-time state space model denoted by 

[A -1 B C D G k] 

Using the relationship, B = -A~ l B, K = -A~ l K, the forward-time system and ob- 
server/Kalman filter model represented by 

[A B C D G K] 
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can be easily computed. Note that G is a backward-time observer gain for the system matrix 
inverse A~ l ,K is the Kalman filter gain, and the matrices A,B,C,D form the state space 
model of the system being identified. 

Experimental Results 

The developed method is applied to data obtained from the truss structure shown in 
Fig. 1 at NASA Langley Research Center. The L-shaped aluminum truss is oriented such 
that its long section is in a vertical direction extending 90 inches. The short section, 20 
inches long, is horizontal and is clamped at the free end to a steel rigid plate. The square 
cross section is 10 inches by 10 inches. Two cold air jet thrusters, located at the beam tip, 
serve as actuators for excitation and control. Each thruster has a maximum thrust of 2.2 
lbs. Two servo accelerometers located at a corner of the square cross section provides the 
in-plane tip acceleration measurements. In addition, an offset weight of 30 pounds is added 
to enhance the dynamic coupling between the two principal axes, and to lower the structure 
fundamental frequency. For identification, the truss is excited using random inputs to both 
thrusters. The acceleration signals are filtered by a three-pole analog Bessel filter with a 
break frequency of 20 Hz. The input-output signals are sampled at 250 Hz and recorded for 
system identification. The following examples illustrate the application of system identifi- 
cation method developed in this paper using both backward-time ARX and backward-time 
ARM AX model approaches. All results are obtained using the System/Observer/Controller 
Identification Toolbox 14 developed at NASA Langley Research Center. 

Identification Using Backward-time ARX Model 

An input-output data record of 2000 points is used in the computation. The order of the 
ARX model is set to be 10, i.e., p = 10. Equation (15) is used to compute the backward- 
time observer Markov parameters from which the Markov parameters of the backward-time 
system are computed. Application of ERA produces a backward-time state space model 
which is then converted to a forward-time representation. Recall that in the backward-time 
model, stable system modes appear unstable, whereas the noise modes remain stable. When 
converted to a forward-time representation, the system modes become stable and the noise 
modes will become unstable and can be eliminated. In other words, this process is simply 
model reduction by modal truncation where the unstable noise modes (in the forward-time 
representation) are discarded. The final reduced order model so obtained contains only the 
stable system modes that are identified by the method. 
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Figures 2a and 2b show the reconstructed responses for both outputs using a 6-order 
state space model obtained by this approach. In these and the following figures, the solid 
curves are the actual responses that are used in the identification procedure, and the dashed 
curves represent the reconstruction. The quality of the reconstruction shows how well the 
identified model represents the actual system in terms of its ability to reproduce the actual 
responses when driven by the same input time histories. For convenience, the backward-time 
identification method presented in this paper is referred to as backward-time OKID, whereas 
the method formulated in Ref. 10 is referred to as forward-time OKID. When appropriate, 
the type of model used (ARX or ARMAX model) is also indicated. As a baseline comparison, 
the reconstruction plots using the 18-th order model obtained by the forward OKID method 
with the same set of data are also shown in Figs. 3a and 3b. 

Identification Using A Backward-time ARMAX Model 

Recall that the primary purpose of using the ARMAX model structure is to improve 
results obtained from the ARX model alone by modelling the noise dynamics. In order 
to show the relative merit of the iterative procedure using the backward-time ARMAX 
approach, the data record that was used to obtain the results presented in Figs. 2 and 3 
is cut short to 1000 points, resulting in a worsened model obtained by the ARX approach. 
This is reflected in its poor ability to reproduce the data as shown in Figures 4a and 4b. 

The iterative procedure is then applied to model the noise dynamics and correct the 
original model. Figures 5a and 5b show the reconstruction plots after one cycle of iteration. 
The results deteriorate but they are explainable since the procedure starts with the assump- 
tion that the least-squares residual obtained with the original ARX solution is the same as 
the whitened residual of the ARMAX model. Figures 6a and 6b show the reconstruction 
plots after two cycles of iteration. The identification results are clearly improved. Further 
improvement can be obtained with additional iterating cycles. However, experiences have 
shown that later cycles produce little improvement in the results. 

As another comparison of the relative performance between the different methods, Ta- 
ble 1 shows the identified frequencies and damping factors for each of the respective cases 
considered above. In theory, if one uses a backward-time approach then the unstable modes 
in the forward-time representation can be attributed to noise and eliminated from the final 
models. Comparison with forward-time model results, however, suggests that these modes 
could in fact be actual modes, but they are not identified by the backward method due to 
factors such as slight non-linearities, non-ideal noise conditions in the test data. Inclusion of 
these modes in the model would bring the final state space model order from 6 to 18, which 
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is the model used to generate Figs. 2a and 2b. Table 1 also shows the orders of the state 
space models that are used to compute the reconstruction plots shown previously, and the 
norm of the residual using the identified model for each respective method. On the other 
hand, a careful examination of the data shows that there is little improvement in the quality 
of the reconstruction by including these modes in the final state space models. Taking into 
account that the data record is short, this is a case where it is not clear to distinguish a 
system mode which is weakly present in the data from a noise mode. The backward-time 
approach, however, gives a strong indication about the true system modes. This is evident 
from the consistencies in the modes identified from both backward-time and forward-time 
approaches. 


Concluding Remarks 

A system identification procedure using a backward-time approach has been formulated. 
A fundamental contribution of this paper is the connection of the state space model to 
backward-time finite difference models, which include both a backward-time deterministic 
auto-regressive model and a backward-time stochastic autoregressive moving average model. 
In the former case, the connection is established in terms of a backward-time observer, and 
in the latter case it is established in terms of a backward-time observer and a Kalman fil- 
ter. The deterministic formulation is considerably simpler in that it is a linear parameter 
estimation problem. The stochastic approach, which requires modelling of the noise dy- 
namics, can provide improved results over those obtained from the deterministic approach. 
This is accomplished at the expense of an iterating solution to a non-linear parameter es- 
timation problem. A key feature of the backward-time models is that stable system modes 
appear unstable whereas the noise modes remain stable. This offers an opportunity to dis- 
tinguish system modes from noise modes. Since the final state space models contain only 
the stronger stable modes, the backward method generally yields models of smaller order 
when compared with models obtained from a forward method. It is advantageous to use the 
backward method in conjunction with a forward method, so that the effective modes in the 
final model can be better evaluated. Since the model structure of the backward-time model 
is very different from that of the forward-time model, examination of the results provided 
by the two approaches can produce valuable insights into the system being identified. 
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Table 1: Identified frequencies and damping values of the truss structure. 



Mode 

OKIDJt 

OKID-b* 

OKID-b* 

OKID_b* 



ARX 

ARX 

ARX 

ARMAX 

N= 


2000 pts 

2000 pts 

1000 pts 

1000 pts 

Freq. 

1 

5.87 

5.89 

5.92 

5.89 

(Hz.) 

2 

7.29 

7.30 

7.33 

7.33 

3 

19.91 

18.82 

18.54 

15.86 


4 

46.60 

46.64 

46.26 

46.14 


5 

48.61 

48.57 

48.61 

48.64 


6 

64.41 

64.54 

64.15 

60.98 


7 

73.76 

73.76 

73.43 

74.37 


8 

105.63 

105.51 

104.01 

99.83 


9 

110.90 

110.97 

111.43 

113.52 


Damp. 

i 

0.65 

0.61 

0.64 

0.50 

(%) 

2 

0.47 

0.29 

-0.17 

0.89 

3 

61.52 

0.91 

0.76 

21.73 


4 

2.08 

-2.32 

-2.93 

-0.85 


5 

1.53 

-0.77 

-0.95 

-0.63 


6 

9.33 

-0.944 

-7.80 

-0.56 


7 

0.74 

-1.61 

-2.19 

-0.55 


8 

10.24 

-10.21 

-9.29 

-11.27 


9 

5.12 

-5.26 

-5.40 

-2.11 


State Space 
Model Order 

18 

6 

4 

6 

Residual 

0.256 

0.317 

0.673 

0.318 


f Forward-time OKID, * Backward-time OKID 
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Fig. 1: Truss Structure Test Configuration. 
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Time (sec.) 

Fig. 2a: Actual vs. reconstructed response for the first 
output ( — Test, - - Reconstruction): Backward-time 
OKID via an ARX model of order p = 10 with 2000 
data points. 



Time (sec.) 

Fig. 2b: Actual vs. reconstructed response for the sec- 
ond output ( — Test, - - Reconstruction): Backward- 
time OKID via an ARX model of order p = 10 with 
2000 data points. 
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Fig. 3a: Actual vs. reconstructed response for the first out- 
put ( — Test, - - Reconstruction): Forward OKID via an 
ARX model of order p = 10 with 2000 data points. 



Time (sec.) 

Fig. 3b: Actual vs. reconstructed response for the second 
output ( — Test, - - Reconstruction): Forward OKID via 
an ARX model of order p = 10 with 2000 data points. 
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Fig. 4a: Actual vs. reconstructed response for the first out- 
put ( — Test, - - Reconstruction): Backward OKID via 
an ARX model of order p = 10 with 1000 data points. 



Fig. 4b: Actual vs. reconstructed response for the second 
output ( — Test, - - Reconstruction): Backward OKID 
via an ARX model of order p = 10 with 1000 data points. 
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Time (sec.) 

Fig. 5a: Actual vs. reconstructed response for the first 
output ( — Test, - - Reconstruction): Backward OKID 
via an ARMAX model of order p = 10 with 1000 data 
points after one cycle of iteration. 



Fig. 5b: Actual vs. reconstructed response for the second 
output ( — Test, - - Reconstruction): Backward OKID 
via an ARMAX model of order p = 10 with 1000 data 
points after one cycle of iteration. 
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Time (sec.) 

Fig. 6a: Actual vs. reconstructed response for the first 
output ( — Test, - - Reconstruction): Backward OKID 
via an ARMAX model of order p = 10 with 1000 data 
points after two cycles of iterations. 



Time (sec.) 

Fig. 6b: Actual vs. reconstructed response for the second 
output ( — Test, - - Reconstruction): Backward OKID 
via an ARMAX model of order p = 10 with 1000 data 
points after two cycles of iterations. 
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